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Direct simulation of polymer drag reduction 
in free shear flows and vortex dipoles 

By P. Orlandi 1 G.M. Homsy 2 and J. Azaiez 2 


One of the most efficient techniques for drag reduction is the injection of polymers 
near a wall which can achieve a reduction in drag up to 80%. Several experimental 
observations tend to indicate that polymers modify the turbulence structures within 
the buffer layer and show that the changes consist of a weakening of the strength of 
the streamwise vortices. In this paper, we investigate the effects of viscoelasticity on 
two different types of flows: the vortex dipole impinging walls to model streamwise 
vortices in a turbulent boundary layer and the mixing layer that represents free 
shear flows. For this purpose, we examined three different rheological models: the 
Oldroyd-B model, the Jeffrey’s corotational model, and the FENE-P model. 


1 Introduction 

Evidence of drag reduction both by passive and active control has been observed 
experimentally, but a clear explanation of the mechanisms responsible for this re- 
duction has not been given, mainly because the experimental observations cannot 
describe all the details of the flow. In the literature, there is a large number of pa- 
pers devoted to the experimental study of the wall layer structures in drag reducing 
flows. This literature is summarized in the review article of Tiederman (1989). The 
main conclusion that we can draw from these papers is that drag reduction is due 
to modifications of the wall layer structures, particularly in the buffer region, the 
most active region in wall bounded flows. From flow visualizations, Oldaker and 
Tiederman (1977) concluded that the polymer solution inhibits the formation of low 
speed streaks and that, when these are formed, their spacing in wall coordinates 
increases with polymer concentration. 

In the last decade, direct simulations of wail bounded flows have been a very 
useful tool for a deeper understanding of turbulence structures and their role in 
wall friction. Free shear flows of viscoelastic fluids on the other hand, have not 
received as much attention as bounded flows did, and to our knowledge, only few 
and inconclusive experiments have been conducted. 

Due to the lack of a universal constitutive equation that describes most common 
viscoelastic behaviors, numerical studies of non-Newtonian fluids have been limited 
to special types of polymeric solutions in simple flows. In spite of these limitations, 
many results showing important effects of viscoelasticity for different types of flows 
have been reported (Tiederman(1989)). In a previous paper (Orlandi 1991), a 
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heuristic relation between the polymer stresses and the flow properties in a channel 
flow was used to investigate how the flow structures change in a situation where 
drag reduction was achieved. The model was not tensorial invariant and was not 
related to the molecular structure of the polymers. 

In this paper, we wish to initiate a systematic study of the effects of polymers 
on the flow structures using rheological models based on transport equations for 
the stresses of the polymers. In particular, we studied the case of vortex dipole 
in the presence of walls. A vortex dipole models the streamwise vortices in the 
buffer region of a turbulent channel, which are responsible for turbulence produc- 
tion and turbulent drag. This was shown numerically in a quasi two-dimensional 
simulation by Orlandi and Jimenez (1991). We have also examined the case of the 
roll-up of a two-dimensional mixing layer. Our interest in the mixing layer problem 
was motivated by the results of linear stability analysis reported by Azaiez and 
Homsy (1992). This analysis shows that, in a special limit of the elasticity number, 
viscoelasticity reduces the instability of the flow. 

We considered three viscoelastic models: the Jeffrey’s corotational model, the 
Oldroyd-B model, and the FENE-P model which describes the rheological proper- 
ties of dilute polymeric solutions. The flows governed by these rheological models 
are characterized by three dimensionless numbers: The Reynolds number, Re, the 
Weissenberg number, We, a dimensionless measure of the elasticity of the fluid and 

the coefficient A" = — — , a ratio of the solvent viscosity t?c and the polymeric 
VS + VP * J 

contribution to the shear viscosity, rjp. In the case of the FENE-P model, a third 

parameter, 6, related to the nature of the spring used to model the macromolecule, 
is introduced. 

The first interesting physical aspect that comes out of the study of these three 
models is that the formation of small scales is enhanced and that these small scales 
are rapidly dissipated. This unexpected phenomenon needs an experimental confir- 
mation, and we hope that in the near future two-dimensional experiments will be 
performed in order to both understand in details the behavior of dilute polymeric 
solutions and validate the results of our two-dimensional viscoelastic simulations. 

While for the mixing layer we obtained results in a quite large range of the 
Weissenberg number, only results at small Weissenberg numbers could be obtained 
for dipoles impinging walls because of the special flow topology: with the FENE-P 
model, the polymeric stresses reached very large values at the stagnation point, and 
the calculations were diverging for We of order 1. At smaller values of We, the 
calculations did not differ from the Newtonian case because the flow stresses were 
able to overcome any possible polymeric effect. 

2 Physical and numerical model 

We used a vorticity-streamfunction formulation for the Cauchy’s momentum 
equation. This equation is closed through evolution equations relating the stress 
tensor to the shear rate tensor. In all the subsequent analysis, the stress tensor is 
written as the sum of two stresses. 
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T ij =T } (Ky ij +(l-K)a >j ) (1) 

We fixed K = 0.5. 

The first term in Eq.(l) reflects the contribution of the Newtonian solvent and 
is proportional to the shear rate tensor 7 , the second one represents the polymeric 
contribution and is proportional to the tensor a. The tensor a satisfies different 
evolution equations corresponding to the various rheological models we are exam- 
ining. 

Eq.(l) leads to the following vorticity equation 


du k„ 1 — fc /da \2 dan da u - a 22 N ^ 

^ = -j(*,vo + -v u + — (-aif - &r - < 2 > 

where the nonlinear terms have been expressed as a Jacobian. The streamfunction 
is evaluated from the Poisson’s equation = ~ u> - 

In the case of the mixing layer, we imposed periodic conditions in the streamwise 
direction. In the transverse direction, free slip conditions are imposed in the finite 
difference code while periodic conditions are imposed for the spectral code. In the 
case of the spectral code, the domain was made large enough in the transverse 
direction to ensure that all functions go to zero at the periphery of the domain. 
Non-slip conditions are imposed at the wall for the vortex dipole. 

The stresses are governed by a set of partial differential equations which constitute 
the rheological model. We have considered several models: the first one is the 
Jeffrey’s corotational model which is obtained by formulating the equations of state 
in a frame translating with the fluid and rotating with the local angular velocity of 
the fluid. The equation describing this model is: 

a ij + We^=j ij (3) 


da\\ — 022 


where We — — is the Weissenberg number. 

0 

Da t j Dciij /j\ 

- - — — =rr + + Ujkdik (4) 

Dt Dt 

is the Jaumann derivative. 

The rheological equation for the Oldroyd-B model is: 

a tj + We 6 -^- = 7i> (5) 

ScL 

The upper convective derivative is related to the Jaumann derivative through 

ot 

the relation: 


+ iikatj) 


( 6 ) 
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FIGURE 1. Time evolution of the polymeric stresses at We = 5. Oldroyd; 

. . . Jeffreys; FENE-P 


To the contrary of what is observed for the Jeffrey’s model, the Oldroyd-B model 
has a very large dependence on the Weissenberg number as will be discussed in 
the next paragraph. We found that at certain values of We, the stresses grow 
indefinitely in time and the growth rate is flow dependent. 

The Oldroyd-B model gives a steady state elongational viscosity that goes to 
infinity at a finite elongational rate. This unlikely behavior results because the 
Hookean dumbbell model permits infinite extension. In order to avoid this unreal- 
istic behavior, a Warner law is used instead of the Hook law leading to the FENE-P 
model. This model is described by Mackay and Petrie (1989) and the rheological 
equation is: 

a 0 Z + We 6 -^- = 7 0 + ^(1 + a, } We) (7) 


Yl W C 

where Z = 1 + -r (H a,,), n = 2 for a two dimensional flow. In order to 

b n 7 

facilitate the numerical solution of Eq.(7), - ^ was derived from the transport 


equation for 


Dt 


The set of governing equations have been discretized by a finite difference scheme 
second order accurate in space and in time. The main characteristics of the numer- 
ical method is that it discretizes by the Arakawa scheme which conserves 

energy, enstrophy, and skew symmetry in the inviscid limit. The polymer contribu- 
tion to the stresses have been located at different grid positions with an and a 22 at 
the center of the cell and at the same location of u> and 0. This choice allows a 
very compact formulation of the right hand side of Eq. (2) and leads to the solution 
without imposing boundary conditions for an and a 2 2 . On the other hand, a J2 has 
been taken equal to zero both on free-slip and no-slip walls. The stream function 
has been solved directly by using the Fourier transform in the periodic direction. 
Numerical simulations have been performed on the CRAY YMP84 at NASA Ames 
and the CPU required for 1000 time steps on a 128 x 128 grid was 150 seconds. 
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FIGURE 2. Time evolution of the polymeric stresses at b = 5. We — 5; 

We = 20; ...We = 50 


3 Results 


1 Mixing layer 

The mixing layer simulations have been performed in order to study the effects of 
viscoelasticity on the vorticity field. Numerical simulation using a pseudospectral 
method were performed in order to compare with the finite difference results. The 
good comparison between finite difference and pseudospectral simulations will be 
presented elsewhere. 

The initial base state is given by a hyperbolic tangent velocity profile on which 
we superposed a cos perturbation on the stream function with the maximum at the 
centerline y = 0. The most unstable wave number, a — 0.44, has been used in all 
the simulations. The extension of the domain in the streamwise direction is set to 
—6 with 8 the momentum thickness of the mixing layer. The vertical extension has 

^ 8uq 

been heuristically fixed equal to 86. The Reynolds number is defined as Re = — — 

where u 0 = (U\ - J7 2 )/2- U\ and U 2 are the free-stream velocities. 

The growth of the maximum value of( an — 0 , 22 ) an d of a \2 obtained for the three 
viscoelastic models axe shown in Fig. 1. At We = 5, the stresses are smaller for the 
Jeffrey’s corotational model than for the other two models; as a consequence, the 
vorticity field does not differ appreciably from the Newtonian case. In the case of 
the Oldroyd-B model, the stresses grew rapidly due to the unrealistic characteristic 
of this model that allows the macromolecules to extend infinitely. 

We suspected that the divergence of the codes for the Oldroyd-B model at high 
We to result from some numerical instability. By refining the mesh and varying 
the time step for both the finite difference and the pseudospectral codes, we found 
that, at the same values of We, the solution was always diverging at the same point 
in time. We have shown that the divergence was caused by the change of sign of 
the eigenvalues of the system of nonlinear equations. The change of sign occurred 
at certain values of We, and this value is flow dependent. 

In the case of the FENE-P model, the dependence of the flow on the Weissenberg 
number and on the parameter b has been studied by performing different simulations 
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FIGURE 3. Time evolution of the polymeric stresses at We = 20. 6=1; 

6 = 5; . ..6= 10 

for the mixing layer at Re = 100. Fig. 2 shows how the polymeric contribution to 
the stresses depend on We at 6 = 5. The simulations at We = 5 and We — 20 
were performed on a 128 x 128 grid. With this grid, at We — 50, the high values of 
the stresses produced wiggles in the vorticity field at t « 40; these wiggles partially 
persisted with a grid twice finer and disappeared with a 384 x 384 grid. Fig. 2 also 
shows that the effect of the Weissenberg number on a \2 are large at later times and 
that large variations of the maximum value of( an — 022 ) occur at low rather than 
at high We values. Fig. 2 shows also that by increasing the Weissenberg number 
from We = 20 to We = 50, the variations are less important than those obtained 
by taking the We from 1 to 20. 

Fig. 3 shows the dependence on the parameter 6. We notice that the maxima are 
more dependent on 6 than on We and that the maximum values are higher than 
those observed in Fig. 2. This behavior is consistent with the fact that for very 
large values of 6 the Oldroyd-B model is recovered. 

A clear picture of the evolution of the flow is obtained from the analysis of the 
vorticity and of the stress field . Vorticity contour plots at three different times 
and for three different We values are presented in Fig. 4. We remark that at low 
We, the flow does not differ from the Newtonian case while at higher We, more 
noticeable changes appear in the structure of the mixing layer. Very interesting is 
the occurrence of intense gradients in the braids. These gradients are convected in 
the roll-up region and persist for a longer time, producing a faster and more intense 
roll-up as shown at We = 50. 

an — a 22 is the other quantity contributing to the modification of the vorticity 
field and can be measured experimentally. Fig. 5 shows a large increase of its 
maximum value when We goes from We = 1 to We = 20. On the other hand, 
when the Weissenberg number increases from We = 20 to We — 50, the value 
does not change significantly, then we can expect minor changes in the vorticity 
field. From Fig. 4 and Fig. 5 we can draw a first conclusion that the effect of the 
Weissenberg number is to concentrate the modifications by polymers at small scales. 
From contour plots of vorticity and an — 022 , not reported in this paper, we have 
observed that by increasing 6, high peak values are produced and are localized in 
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t = 15 t = 30 t = 45 

FIGURE 4. Vorticity contour plots. b=5. 

wider regions, and thus weaker effects on the mixing layer roll-up are observed. 

At this stage of the work and without having any experimental observations with 
which to compare our results, we do not wish to go into a deeper analysis of the flow 
to understand the mechanisms responsible for this rapid transfer of energy from the 
large scales to the smaller scales. We think that the results we obtained will be of 
more interest if they can be checked through experimental studies. The comparison 
between numerical and experimental predictions will help in understanding the 
physics behind the effects of viscoelasticity on the structure of the flow. 

2 Vortex dipole impinging walls 

This two-dimensional flow has been considered because it represents the stream- 
wise vortices in a turbulent boundary layer, which are responsible of turbulent 
drag and turbulence production. In the present two-dimensional simulation, we 
have introduced a Lamb dipole rather than a single vortical structure because the 
dipole moves towards the wall by self-induction while a single vortex must be ad- 
vected towards the wall by an external field. As a first case, we have considered 
the interaction with a free-slip wall to investigate the distributions of the polymeric 
contribution to the stresses generated at the moment of the impact. We expect that, 
at the stagnation point, a large increase of the polymeric stresses is obtained since 
this is a point where strong elongations of the polymers molecule occur. This large 
increase of the polymers stresses is the cause of the difficulty to perform simulations 
at high Weissenberg numbers. 
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t = 15 t = 30 t = 45 


FIGURE 5. an — <222 contour plots. b=5. 

The simulations were performed on a 192 x 192 grid points in a domain extending 
in the horizontal and in the vertical directions 4 dipoles radii. Solutions with the 

FENE-P model were obtained only for We < 3 and b < 1 at Re = = 50 (Ud 

v 

is the translation velocity of the dipole, a is the dipole radius). The simulations 
were performed both with free-slip and no-slip walls. In the case of the dipole 
impinging a free-slip wall, the vorticity field, at the moment of the impact, is shown 
in Fig. 6. Similarly to what has been observed in the mixing layer, this figure 
shows that an — 022 increases more than a\ 2 ] the maximum values are reached in 
the region of maximum normal stress and low vorticity. Although large polymeric 
stresses are obtained, the dipole does not change their shape appreciably, and this 
could be a consequence of the fact that dipolar structures are very stable to large 
perturbations. 

When the dipole impinges a non-slip wall, Fig. 7 shows that the vorticity field 
does not depend strongly on the Weissenberg number. There is, however, a slight 
decrease of the magnitude of the vorticity level as compared to the Newtonian case, 
but this decrease is not sufficient to change the type of vortex rebound. We did not 
perform a systematic study of the We effects. 

4 Conclusion 

This preliminary study testing different viscoelastic models allowed us to predict 
the effects of the introduction of polymers on the flow structure. By comparing 
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(a) (b) (c) 

FIGURE 6. Contour plots of (a) vorticity, (b) an — a^ and (c) a \2 for dipole 
impinging a free-slip wall at Re — 50 , We = 2 and b = A. 



FIGURE 7. Vorticity Contour plots (a) newtonian, (b) We = .5,6 — 1. and (c) 
We = 3., b = 1. for dipole impinging a non-slip wall at Re — 50. 


the solutions obtained by the finite difference scheme and those obtained by a 
pseudospectral method, we have shown that the solutions did not diverge in time 
for numerical reasons. Two-dimensional simulations revealed that the cause of the 
exponential growth in the Oldroyd-B model at certain W e is caused by a change in 
the sign of the eigenvalues of the an stress equation. 

Using the FENE-P model, we were able to conduct numerical simulations at rea- 
sonably high Weissenberg numbers. In the case of the mixing layer, we noticed some 
changes in the structures of the roll-up. This is an interesting phenomenon since, 
in a space developing mixing layer just after the splitter plate, the initial instability 
could influence the growth of the mixing layer. The present simulations show the 
large potential in using polymers for the control and possibly the modification of 
the structures of the mixing layer. We hope to pursue this study by examining the 
effects of polymers on vortex stretching and streamwise vortices in the mixing layer. 
In the case of the vortex dipole, the changes on the vorticity field are negligible. 

We believe that the results reported in this study should be pursued by more 
detailed simulations and by experimental studies that will help in understanding 
the mechanisms responsible for the changes observed in the case of the mixing layer 
with the FENE-P model. We hope that the present study will stimulate more 
interest in non- Newtonian flows among both theoreticians and experimentalists. 
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